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Q ' boundaries is considered. The stationary state, which is known to be constructed in 

a matrix product form, is studied by applying the theory of g-orthogonal polynomials. 
Using a formula of the g-Hermite polynomials, the average density profile is computed 

"^ I in the thermodynamic limit. The phase diagram for the correlation length, which 

r^^ ■ was conjectured in ||l|(J. Phys. A 32 (1999) 7109), is confirmed. 
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t^ 1 Introduction 



The one-dimensional partially asymmetric simple exclusion process with open 



The one-dimensional asymmetric simple exclusion process (ASEP) ^ ^ is a system of 



^ \ particles which hop preferentially in one direction on a one-dimensional lattice with hard- 

core exclusion interaction. The ASEP has been studied extensively since it is one of the 
few models which show rich non-equilibrium behaviors and is exactly solvable [Q. Besides, 
the ASEP has applications to many interesting problems such as the hopping conductivity, 
growth processes and the traffic flows 0. 

5^ \ In this article, we consider the stationary state of the ASEP with open boundary 

conditions. That is, the system is connected to particle resevoirs at boundaries. The case 
where particles can hop only in one direction, which we refer to as the "totally asymmetric" 
case in the sequel, was solved in |§|, • The current and the density profile were calculated 
exactly in the thermodynamic limit. The phase diagram for the current and the correlation 
length were identified. The system exhibits phase transitions depending on the parameters 
at the boundaries. Recently the obtained phase diagram was discussed from the point of 
view of the domain wall dynamics |^ . 
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The partially asymmetric case with the open boundary conditions was partially solved 
in IQ. The current was evaluated in the thermodynamic limit. The phase diagram for 
the current was identified. It turned out to be the same as the one obtained by mean- 
field approximation [Q or by employing a plausible assumption [0. The phase diagram 
for the correlaiton length was also obtained by assuming that the correlation length is 
given by the logarithm of the ratio of the largest and the second largest eigenvalue of a 
certain matrix which plays a similar role as a transfer matrix does in equilibrium statistical 
mechanical models. It was shown that the phase diagram has a richer structure than that 
for the totally asymmetric case. The average density profile was, however, not calculated 
in [m. In this sence the obtained phase diagram for the correlation length has remained 
a conjecture. The purpose of this paper is the confirmation of this phase diagram. By 
using the explicit formula for the Poisson kernel of the g-Hermite polynomials, the average 
density profile in the thermodynamic limit is calculated for the partially asymmetric case. 
It turns out that the phase diagram was correctly predicted in |]I|]. 

In this article, we only consider the case where hoppings of particles at the boundaries 
and those at the bulk part of the system are compatible. In other words, when we allow 
the particle input at the left boundary and the particle output at the right boundary, the 
hopping rate to the right is assumed to be larger than that to the left. When hoppings 
at the boundaries and those at bulk is imcompatible, the current becomes zero in the 
thermodynamic limit. The situation seems to be similar to the closed boundary condition 



where particles can not enter or go out of the system [|Tl|]. Off course, when we consider the 
finite chain, the current remains to be positive. We remark that the asymptotic current 
for this case was evaluated in lO . 



The paper is organized as follows. In the next section, the definition of the model is 
given in terms of the master equation. The so-called matrix product ansatz, which gives 
the stationary state in the form of matrix product, is also explained. Some properties of 
the g-Hermite polynomials and the relationship to the matrix product ansatz are explained 
in section ^. The section ^ is the main section of this article. First, the one-point funciton 
is represented in the form of double integrals. Second, the average density profile in the 
thermodynamic limit are summerlized whereas the evaluation of the integrals are relegated 
to Appendices. The phase diagram for the correlation length is identified. The concluding 
remarks are given in the last section. 

2 Definition of Model and Matrix Product Ansatz 

The one-dimensional asymmetric simple exclusion process (ASEP) is defined as follows. 
During the infinitesimal time interval dt, each particle jumps to the right nearest neigh- 
boring site with probability pndt and to the left nearest neighboring site with probability 
PL^t. If the chosen site is already occupied, the particle does not move due to the exclusion 
rule. More than one particle can not be on the same site. Each site can be either empty 
or occupied. The case where particles can hop only in one direction, i.e., the case where 
either pl = or pr = is called the "totally asymmetric" case. The pr = pi case is 



called the "symmetric" case whereas the case where particles hop in both directions with 
different rates will be referred to as the "partially asymmetric" case. In addition, we allow 
the particle input at the left end of the chain with rate a and allow the particle output at 
the right end of the chain with rate /3 (Fig. 1). Here the length of the chain is denoted 
by L. In this article, we restrict our attention to the partially asymmetric case since the 
totally asymmetric case and the symmetric case was already solved in p, ^ and in ||13| 
respectively. The restrictions on the parameters are ^ <Pl < Pr and a, /9 > 0. 

More formally, the process is defined in terms of the master equation. Each config- 
uration of the system is indicated by {ri, r2, . . . , tl} where Tj [j = 1,2,..., L) denotes 
the particle number at site j. Namely Tj = if the site j is empty whereas Tj = 1 if the 
site j is occupied. Let P(ti, T2, . . . , r^; t) denote the probability that the system has the 
configulation {ti, T2, . . . , r^} at time t. Then the time evolution of the ASEP is described 
by the following master equation. 



— P(ri,r2, ... ,TL;t) 
= a(2ri-l)P(0,r2,... ,TL;t) 

L-l 
+ ^{Tj - Tj+i) [pLPin, T2,... ,0, 1, 

i=i 
-PRPiri,T2,... ,1,0,... ,TL;t)] 
+ /?(l-2ri)P(ri,r2,... ,rL_i,l;t). 



.rL^t) 



(2.1) 



For instance, the master equation for L 



d 
dt 
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P(10;t) 
P(ll;t) 
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(2.2) 



One can confirm himself that the dynamics of the ASEP is correctly encoded in the master 



equation (|2.1| ). 

When time t goes to infinity, the system is expected to reach the stationary state. The 
probability distribution in the stationary state will be denoted as P{ti,T2, . . . ,tl). For 
instance, except for the normalization, the stationary state for L = 2 case is the eigenvector 
of the 4x4 matrix in the right hand side of ( p.2| ) with the eigenvalue zero. Explicitly, it 
reads 



■p(oo)' 

P(01) 
P(10) 

P(ll) 



Const. 



j_ 

PR 



J_ 

cP 
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a/3 ^ a ^ /3 
1 



PL 



(2.3) 



In 1^, it was shown that the probabihty distribution of the ASEP in the stationary 
state for general L can be written in the form of the matrix product as 



P(ri,r2,... ,tl) 



Uw\\{{r,D + {l-r,)E)\V), 



(2.4) 



i=i 



where D and E are square matrices and (W\ and \V) are vectors satisfying following 
relations, 



PRDE-pLED = aD + E), 

a{w\E = aw\. mv) = (:\v). 

Here C, is an arbitrary number. If one defines the matrix C by 

C = D + E, 
the normalization Zl is given by 

zl = mc''\v). 



(2.5a) 
(2.5b) 



(2.6) 



(2.7) 



Here, for the case of L = 2, we check that the state (|2.4|) indeed gives the stationary 
state of the process by using the algebraic relations (|2.5a|) and (|2.5b|) . For L = 2, (|2.4| ) 
reads 



■p(ooy 

P(01) 
P(10) 
Pill) 
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'{W\E^\V)' 
{W\ED\V) 
{W\DE\V) 
{W\D^\V) 



f2.8) 



Three components, P(00), P(01), P(ll) can be calculated by simply using (|2.5b|) . On the 
other hand, one computes P(10) first changing the order of matrices D, E by ( |2.5a| ) and 
then using ( p.5b| ). Hence we get 



■p(ooy 

P(01) 
P(10) 
P(ll) 
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^(ei^ + 1 

PR \al3 a 
/32 



(2.9) 



One can compare this expression with ( p.3|) to see that this expression indeed gives the 
stationary state for L = 2 case. One also sees that the arbitrary parameter ( appears in 
the same way for all components, P(00), P(01), P(10), P(ll). Changing the parameter ( 
only changes the normalization Z2. This is true for general L as well. In this article, the 
proof that the state (|2.4| ) gives the stationary state for general L is not given. See . 



We express several physical quantities in the form of matrix products. The one-point 
function {nj)L is defined as the probability that the site j is occupied. In other words, 
{nj)L is the average density at site j. The two-point function {njnk)L is defined as the 
probability that the sites j and the site k are both occupied. Higher correlation functions 
are defined similarly. In the matrix language, they are computed by 

{n,)L = {W\C^-'DC'^-^V)/Zl, (2.10) 

{n,nk)L = {W\C^''DC''~^~'DC'^-''\V)/Zl, (2.11) 

and so on. The current through the bond between site j and site j ' + 1 is defined by 
Jl — P-r(^j(1 ""^j+i)) ~PL{nj{l —rij^i)). In the steady state, the current is independent 
of j and hence is denoted by Jl ■ It is given by 

'^'-^ {w\c^\v) -^^T' ^ ^ 

Once one finds a representation of these algebraic relations, by using the above formula, one 
can in principle calculate the physical quantities such as the particle current Jl, the one- 
point function {nj)L, the two-point function {njnk)L and the higher correlation functions. 
We note that the process has an obvious particle-hole symmetry. When we look at 
holes instead of particles, they tend to hop to the left with rate p^j and to the right with 
rate pl with hard-core exclusion. In addition, they are injected at right end with rate (3 
and they are removed at the left end with rate a. In other words, the process is invariant 
under the changes, 

particle ^^ hole 

a^p (2.13) 

site number j ^^ site number L — j + 1. 

Due to this symmetry, it is sufficient to obtain the density for the right half of the system. 
The density for the left half of the system is obtained by using the above symmetry as 

{n,)L{a,l3) = 1 - K_,+i)i(/3,a), (2.14) 

where the dependence of {nj)L on the parameters a and (3 are eplicitly indicated. 

Before closing the section, we present some simulation results (Fig. 2 and Fig. 3). 
Figure 2 shows the space-time diagrams for several choices of parameters. It is clear that 
the properties of the system crutially depend on the values of the boundary parameters. 
The differeces become more transparent when we consider the particle current or the the 
average density profile of the stationary state. In principle, the stationary state is achieved 
only in the infinite time limit. However, we see from Fig. 2 that the system practically 
goes into a stationary state after some transient time. Hence if we average the density over 
a long time after the transient time, it would be regarded as the average density profile of 
the stationary state practically. The results are shown in Fig. 3. When a is small and f3 



is large, the bulk density is low. It decays sharply near the right boundary. This is called 
the low density phase. Conversely, when a is large and /3 is small, the bulk density is high. 
It decays sharply near the left boundary. This is called the high density phase. When 
a = /5 is small, the low density region and the high density region coexist (coexistence 
line). Finally, when both a and /5 are large enough, the density takes the value 1/2 at bulk 
and decays slowly near both the boundaries. This is called the maximal current phase. 
Our main tasks in the following are to obtain the average density profiles in Fig. 4 exactly. 

3 Representation of Algebra and g-Hermite Polyno- 
mials 

First we introduce some notations for later convenience. We introduce the g-number, 

W = l-g^ (3.1) 

and the g-shifted factorial, 

(a; q\ = (1 - a)(l - aq){\ - aq^) ■ ■ ■ (1 - ag^-^), (3.2a) 

(a;g)o = l. (3.2b) 

We also define 

oo 

{a](l)oo = \\{l~aq^), (3.3) 

for |g| < 1. Since products of g-shifted factorials appear so often, we use the notations, 

(ai, 02, ■ ■ ■ , Ofc; g)oo = (en; g)oo(a2; q)oo ■ ■ ■ {ak, g)oo, (3.4) 

(ai, 02, ■ ■ ■ , ak] q)n = (en; q)n{a2] q)n- ■■ [ak] q)n- (3.5) 



Next a representation of the algebraic relaiton ( ^.5a] ) and ( |2.5bD is given. If we define 



D = l + d, E = l + e, (3.6) 
we see that these relations become 

de — qed = 1 — g, (3.7a) 

{W\e = a{W\, d\V) = b\V), (3.7b) 



where we put 



<? =Pl/pr, (3.8) 

a=^^, b=^^—, (3.9) 



with a = a/{j)R — Pl), P = P/{pr — Pl)- Since < pr < pi, we have < g < 1. 

In this article, we take the following representation for the matrices d, e and the vectors 
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{1}5 











{2}i 











{3} 



{1}^ 



















{3}^ 










(3.10a) 



/ 1 \ 

b 



(iy| = /€e(a| = «: 1, 



\/(?;9)i' \/{(i\q)2 



\V) = K\h), = K 



\/{q\q)i 
b^ 



V '■■ J 

(3.10b) 

The constant k is takes as k^ = {ab;q)oo so that (W^|V^) = 1. It should be noticed that 
there exists another useful representaion of the algebraic relations ( p.5a| ) and ( |2.5b| ). It was 
first given in and was used to obtain the phase diagram of the correlation length in [Q . 
The advantage of the represetation ( |3.10a| ) and ( |3.10b| ) is that the commutation relation 
of the matrices d and e turns out to be a simple diagonal matrix. We have 



de — ed = {1 — q) 
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g 

g2 

g3 



(3.11) 



This fact will play an important role for the calculation of the average density profile in 
the next section. 

Next we list some properties of the continuous g-Hermite polynomials . The proofs can 



be found for instance in ||1J, |15|. The continuous g-Hermite polynomials {if„(x|g)|n 
0, 1, 2, ... } are defined by the three term recurrence relation, 



Hn+i{x; g) + (1 - g")if„_i(2;; g) = 2xif„(x; g), 
with the initial condition, 

/7„i(a;;g) = 0, H^{x-q) = l. 

They are explicitly given by the formula. 



H^{cose\q) = Y, 



(?;?) 



fc=0 



iq;q)kiq;q)n-k 



" i{n-2k)e 



(3.12) 
(3.13) 

(3.14) 



The orthogonality relation reads 

/ Hn{cos9\q)H^{cose\q){e^'\e-^''-q)^d9 = I-k 
Jo 

The generating function is also known and is given by 

Y^ Hn{cos9\q) ^^ _ 1 






n=0 



(g;?). 



(Ae*^ Xe-^O; q). 



(3.15) 



(3.16) 



for |A| < 1. To calculate the average density profile, we also need the so-called Poisson 
kernel, 



'-1 

oo 



(r2;g), 



E Hn{cos6\q)Hn{cosLp\q)r'^ _ 



n=0 



(?;?)» 



(3.17) 



This formula is called g-Mehler formula in the mathematics literature. 
Here we notince that, if we introduce 

Pn{x) = Hn{x\q)/^/{q;q)n, 
the three-term recurrence relation is rewritten into the form, 

{1}^ 

{1}^ {2}i 
{2}^ {3}^ 



(3.18) 





Po{x) 




Po{x) 




Pi{x) 

P2ix) 


= 2x 


Pi{x) 

P2ix) 











(3.19) 



In other words, \p{x)) = ^ {po{x) , pi{x) , . . . ) is formally an eigenvector of the matrix d + e 
with eigenvalue 2x. This is the basic relationship between the representaion of the algebra 
( p.5a| ),( p.5b| ) and the theory of g-orthogonal polynomials. Finally, the completeness of the 
continuous g-Hermite polynomials reads 



27r 



d9{e" 



q)oc>\p{cos9)){p{cosi 



(3.20) 



4 Calculation of Density Profile 



In this section, the average density profile is calculated by using the formula (|3.17|) of the 
continuous g-Hermite polynomials. We first recall the asymptotic behaviors of the normal- 
ization Zl and the current in the thermodynamic limit J = lim^^^oo Jl- The asymptotic 
expressions for Zl were given in |^ and are summarized as follows: 

• For phase A (low-density phase; a > 1 and a > b ; a < ^ and a < (3) 



Zr 



{b/a; q). 



-[(l + a)(l + a-^)]^ 



(4.1) 



For phase B (high-density phase; 6 > 1 and a < 6 ; /3 < | and a > /?) 

-[(l + 6)(l + rl)]^ (4.2) 



{a/b; q). 

For phase C (maximal current phase; 0<a,6<l;/5>| and a > |) 

^ _ (a6;g)oo(g;g)L4^+' 



(a -a i)(a ^■,q)^L^^^ ^ ^^^^ ^ „-imL-i 



(4.3) 



On the coexistense hne (a = 6>l;a = /5<|) 

■[(l + a)(l + a-^)r\ (4.4) 



(g;g)oo 

Using (|2.12|) , the current in the thermodynamic hmit is readily computed as 

• For phase A (a < | and P > a) 

J = {pR-pL)a{l- a), (4.5) 

• For phase B (/3 < ^ and a > [3) 

J = {Pr-Pl)P{1-P), (4.6) 

• For phase C (a > | and /? > |) 

J = ^-^^. (4.7) 
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The phase diagram for the current is depicted in Fig. 4. 

Now we turn to consider the average density profile. As you can see from the figures 
in Fig. 4, the average density is almost constant at bulk part except on the coexistence 
line. Hence we are interested in the average bulk density and how the density decays 
near the boundaries. When the average density decays like e~^f^ with r distance from 
the boundary, we refer to S, as the correlation length in this paper. As for the average 
density near the boundares, it is sufficient to compute the density near the right boundary 
due to the symmetry (p.l3|) . The relation (|2.14| ) enables us to know the average density 



profile near the left boundary. On the other hand, the coexistence line should be treated 
separately. Since it is easier to calculate the density difference than the density itself, we 
rewrite the density at site j as 

L-l 

(%)i = Y^{'<^k)L - K+i)l) + K)l- (4.8) 

k=j 



At the right boundary, we have 






IZ 



L-1 



/3 Zl 
■i (L ^ oo) (4.9) 



Using (|4.5| )-( ^r7|) , the density at the right boundary is easily calculated. 
Next we notice that 

1 



nk)L - {nk+i)L = —{{W\C'-'DC'^-'^\V) - {W\C'DC'^-'^-'\V)) 



Zl 



Zl 



.{W\C^~\de - ed)C^-^-^\V). (4.10) 

Zl 

Now one can represent (W\C^^^{de — ed)C^^^^^\V) in the form of double integrals. The 
calculation proceeds as follows. First we notice that 

{W\C''-\de - ed)C''-''-^\V) 



X c(a|C'=^Xcos^))(p(cos^)|(t/e - ed)C^~^-^\p{co^^)){p{cos^)\h)^ 



= {ab;q)Uq;q)l f f- T ^(e2^^ e-2^^ e^^^ e'^^^; g)oo[2(l + cos^)]'=-^[2(l + cos^)]^-'=-^ 
Jo 27r Jo 27r 

X c{o.\p{coii0)){p{co?,6)\{de — ed)\p{cosip)){p{co?,ip)\b)c- (4-11) 

Here the formula (|3.16| ) gives 

Mpi^ose)) = ^^^.e^^-.e.^^j (p(cos^)|6). = J^^T^^T^- (4-12) 

The remaining term {p{cos9)\{de — ed)\p{cos (p)) can also be computed by using the fact 
( PH ) and the formula (|3T7|) . We see that 

oo 

{p{cos9)\{de — ed)\p{co?,(p)) = (1 — g) / ^Pn(cos6')p„(cosyj)g" 

Hn{cos6\q)Hn{cos'^\q)q'^ 



n=0 

oo 



;i-'z)E 



n=0 (^' ^)' 

(?;?)oo 



(ge*(^+'^), ge-*(^+^), ge*(^-'^), ge-*(^-^); g)<. 
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(4.13) 



Hence we have 



{W\C'"\de - ed)C^~''~'\V) 



3 r d0 r dip (e2*^e-2*^,e2^'^,e-2^^;g)oo[2(l + cos^)]'=-i[2(l + cos(^)]^-'=-i 
-H;g)oo(g;g)ocy^ 2^X 2^ {ae'O, ae-^^ ge^(^+^), ge-*(^+^), qe^('-^\ qe-^('-^\ 6e^^, foe"^^; g)o, 

(4.14) 

for < a, 6 < 1. Summing with respect to k from j to L — 1, we get 

L-l 

^{{nk)L - K+i)l) 

k=j 

{ab;q)^{q;q)l H d9 H dip (e^^^ e-2^^ e^^^, e'^^^; g) 



1 y/oo 



Zl Jo 27r Jq 2% (ae^^ ae-^^ ge*(^+'^), ge-'(^+'^), ge*^^-'^), ge-*(^-'^), fee^"^, fee"^"^; g)c 

[2(1 + cosg)p''^[2(l + cos(^)]^-^' - [2(1 + cosg)]^"^ 
2 cos (/? — 2 cos p 



when < a, 6 < 1. Hence, substitution of this expression into ( fl.8|) leads to the expression 
of the average density profile for < a, 6 < 1. Finally, the analytic continuation of ( [4. 151) 



gives the expression of the density porfile for other values of the parameters as well. This 
is conveniently done by writing (|4.15| ) as a contour integral of the two complex variables 

Zi = e*^ and Z2 = e"^. That is, we rewrite (^4.15|) as 

L-l 

J2{(^k)L - K+i)l) = ^(/i + h), (4.16) 

k=j 

where 

dzi f dz2 



h = -^iab;q)ooiq;q)lo j 



2TTiZi J 27TiZ2 

zi z^\ zi z^^- g)oo[(i + 2i)(i + ^r')]^-' 



{azi, az^ \ qziZ2, qz^ ^z^ \ qziz^ \ qz^ ^z^ \ 62:2, 62:2 ^ q)oo{z2 + z^^ - Zi - z^^Y 
1/7 \S [ d^i [ dz2 

4 J 2'KiZi J 2mz2 



1^1)^1 5^2 5^2 ' 



g)oo[(l + zi)il + z^')y-'[il + Z2)il + z^')Y 



{azi, az^ \ qziZ2, qz^ ^z^ \ qziZ2 \ qz^ ^z^ \ 6^2, bz^ ^ 5)00(2:2 + ^2 ^ - ^1 " ^^i ^) 

(4.17) 

When < a, 6 < 1, the contours of zi and Z2 are both unit circles. For other valures of the 
parameters, the contours are deformed as poles in the integrands move in and out of the 
unit circles. 

The evaluation of the integrals /i and I2 in the thermodjTiamic limit for each phase is 
relegated to Appendices. Basically, the integral Ji gives the density at bulk region whereas 
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the integral I2 gives the density near the right boundary. The results are summerized in 
the following. We thus obtain the phase diagram shown in Fig. 5. Notice that this phase 
diagram was corrrectly predicted in [Q. Therein the correlation length for each phase was 
assumed to be given by the logarithm of the ratio of the largest and the second largest 
eigenvalue of the matrix C. 

• Phase C {a> 1/2 and (3 > 1/2; < a, 6 < 1) 

In this phase, the average density at bulk is 1/2. The average density decays near 
the right bounary as 

K)l = ^--^. (4.18) 

Here and in the following, we set / = L — j + 1. The density decays algebraically and 
hence the correlation length is infinite. This average density profile is exactly the 
same as that for the totally asymmetric case. The density decay at the left boundary 
can be obtained from the symmetry relation ( |2.14D . 

• Phase Ai {a < (3 < a/[{l — a)q + a] and (3 < 1/2; aq < b < a and 6 > 1 ) 

The average density at bulk is a. The density near the right boundary decays expo- 
nentially as 



{nj)L = a - 



(6 ^g, g;g)oo a(l -a) 



(a ^b ^q,ab ^q;q)c 
with the correlation length 



P{l-P) 



2/3), (4.19) 



e-^ln^^. (4.20) 

a[l — a) 

On the other hand, the density near the left boundary takes the constant value a. 
This fact is obtained by conbining the average density profile near the right boundary 
for the high-density phase below and the symmetry ( p.l3|) . 



Phase A2 (g/(l + q) < a < 1/2 and P > 1/2; 1 < a < g"^ and 6 < 1) 

In the bulk region, the average density takes the constant value a. On the other 
hand, the density profile near the right boundary decays exponentially as 

_ , (g - &)(1 - ab){abq, g-^bq; g)oo(g; g)^ [4a(l - a)f 
^""'^^ " {a-mb-maq,a-^q,bq;q)o. ^ll ' ^^ ^ 

with the correlation length 

r^ = -ln4[a(l-a)]. (4.22) 

But the decay is not purely exponential but with algebraic corrections. At the left 
boundary, the density takes the constant value a. 
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Phase A3 { P > a/[(l — a)q + a] and a < g/(l + q); a > q ^ and b < aq) 

The average density at bulk is a. Near the right boundary, the density decays 
exponentially as 



(%) 



ilL 



a 



{I - ah){l - {aq)-^) 
(l-6(ag)-i)(l + ag) 

with the correlation length 



'l + aq){l + {aq)- 
(l + a)(l + a-i^ 



r' = iht^ 



q 



\a 



'l-&)qY 



(4.23) 



(4.24) 



This density profile has no correspondense for the totally asymmetric case. At the 
left boundary, the density takes the constant value a. 

Phase Bi {aq/[{l — a) + qa] < (3 < (3 and a < 1/2; bq < a < b and a > 1) 

The average density at bulk is l—(3. As can be seen from the calculation in Appendix, 
The density takes the constant value near the right boundary. 



(n^) 



j/L 



P 



(4.25) 



Since this phase is related to the phase Ai through the symmetry (|2.13| ), the denisty 
decays exponentially near the left boundary with the correlation length 



C 



a(l — a) 
m— ^, 

Pil-P) 



(4.26) 



Phase B2 (g/(l + g) < /3 < 1/2 and a > 1/2; a <1 and 1 < 6 < g"^) 

This phase is symmetric to phase A2 through the symmetry ( |2.13| ). The average 
density at bulk and near the right boundary is 1 — /?. Near the left boundary, the 
density decays exponentially with the correlation length 



r' = -ln4[;3(l-^)]. 



(4.27) 



Phase -B3 (a < bq and b > q ^) (/5 < aq/[l — a + aq] and /3 < g/(l + g); a < bq and 
b>q-^) 

This phase is symmetric to phase A3 through the symmetry ( [2.131 ). The average 
density at bulk and near the right boundary is 1— /?. The density decays exponentially 
with the correlation length 



r^ = In^ 



[P+{l-P)qY 



(4.28) 



near the left boundary. 
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• Coexistence line (a = /? < 1/2; a = b > 1) 

The average density shows hnear profile at bulk, 

{nj)L = a + {l-2a)^. (4.29) 

This is essentially the same as the result for the totally asymmetric case. 

Before closing the section, we present a simulation result for the correlation length to 
show the differece between the phase A2 and ^3. The simulation was done on j3 = 1 
line. For the totally asymmetric case, the system is in the A2 phase on this line. The the 
correlation length is given by ^a2 = l/ln4[a(l — a)]. On the other hand, for the partially 
asymmetric case, the system is in the A3 phase when a < q/{l + q). The correlation length 
is given by ^^13 = 1/ln - ^- ^ . The differences between these two expressions become 
large especially when q and a are small. Especially, as a — ;> 0, ^^2 go^s to zero whereas 
^^3 goes to — 1/lng. In Fig. 6 the simulation result for the correlation length is shown for 
Pr = 1.0, Pl = 0.9, /3 = 10 which corresponds to g = 0.9, /3 = 1. It is clear that the the 
correlatin lenght approaches the finite value as a ^- and is well described by the formula 
for ^^3. 

5 Concluding Remarks 

In this article, we have computed the average density profile of the partially asymmetric 
simple exclusion process with open boundaries. The calculation has been done for a wide 
rage of parameters satisfying < pl < pr and a > 0,/? > 0. The phase diagram for 
the correlation length has been obtained. It has turned out that the phase diagram was 
correctly predicted in the earlier paper |jl[|. In [||], the phase diagram was obtained by 
assuming that the correlation length is given by the logarithm of the ratio of the largest 
and the second largest eigenvalues of the matrix C. The discussions were only for the 
phases with the exponentially decaying profile. In this article, we have not only confirmed 
this fact but also obtained the asymptotic expressions of the average density profile for all 
phases. 

There are two key facts which allowed us to calculate the average density profile exactly 
in the thermodaynamic limit. One is that the commuation relation of the matrices D, E 
becomes a simple diagonal matrix and the other is the formula ( 3.17|) of the g-Hermite 
polynomials. 

There seems to be many possible applications and generalizations of the analysis of this 
article. First, it is possible to generalize the analysis in this paper to the partially asym- 
metric exclusion process on a ring with a single defect particle p!6| , p!7[ . The corresponding 



totally asymmetric case was already solved in |jT6|. Second, the case where pl > pr is also 
interesting. Although the current was evaluated in |]T^, more exact results are desirable. 
Third, it would be interesting to apply the simialr analysis to the multi-species models 



r^, pO, pri, 122]. Compared to the ASEP, much less is known about these models. 
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Several investigatios are now in progress |^. The results about these will be reported 
elsewhere. 
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Appendix A Evaluation of Integral /i 



In this appendix, the integral I 



1-, 



1 / r \ / n3 /" '^'^1 /" '^'^2 



h = -:ia,b;q)ooiq;q)i 



4 J 27cizi J 2Tiiz2 

{azi, az^^, qziZ2, qz];^ z^^ , qziZ^^ , qz^^z^^,bz2, hz^^] q)ooiz2 + ^^^ - Zi - z^^) ' 



(A.l) 



is evaluated. First, for the case where a,b < 1, both of the contours of zi and Z2 are unit 
circles. We have 

""d^ r d^ (e^^^ e-2^^ e^^^, e-2^^; g)oo 



Jo ^TT Jq 



27r Jq 27r (ae*^ ae''^, qe^^+f\ ge-*(^+v'), ge*(^-^), qe-'^^-f\ be''^, te"*'^; g)oo 

xMi±iH£f)l!^. ,A.2) 

2 cos cp — 2 cos 9 

Second, consider the case where a becomes larger one but b is still smaller than one. We 
assume 

a>aq>aq^>---> aq^^"^ > 1 > 05"'"'+^ > ■ • • . (A.3) 

Then the contour of Z2 is still a unit circle but the contour of Zi has to be modified to 
include all poles at zi = aq^ (/c = 0, 1, . . . ,n*^"^) and to exclude all poles at zi = {aq'')~^ 
{k = 0,1, . . . , n^")). Seperating the contributions from poles at zi = aq'' and Zi = [aq'^y^, 
the integral Ji can be rewritten as 

(a) 



J ^ AO) _ iab;q)oc{q;q)i, sr^ (-)^g 

^1 On 2^ 



I ^ (-)V(^-')/'(«'g'',a-2g-2^g)^[Al'^)]^"i 



2a ^ (g;g)fc(a^g'';g)oo 

X f _^ {zlz^^;q)^ ^^^^^ 

Jco 2iTiz2 {aq''+^Z2, aq''+^Z2^, a-'^q-''z2, a-^q-''Z2^, bz2, bz^^; g)oo 

Here the contour Cq denotes the unit circle. The A^'^ s are defined by 

\i'^ = il + cq'){l + c-'q-'), (A.5) 
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for c = a,b and A; = 0, 1, 2, ... . The intgral can be evaluated explicitly by using the general 
formula, 

dz {z^,z-'^;q)^ 2{abcd;q)oo ,. „. 

(A.6) 



Q2mz {az^az ^,bz,bz ^,cz,cz ^,dz,dz ^;q)Qo {q,ab,ac,ad,bc,bd,cd;q)oo 

The contour C is such that it includes all poles of the type fq^ and excludes all poles of the 
type f^^q^^ with f = a,b,c,d and fc = 0, 1, 2, . . . . The parameters a, b, c, d in this formula 
has nothing to do with the a, b, c, d which appear in the rest of this article. This formula 
plays a crutial role in proving the orthogonaliry relation of the Askey- Wilson polynomials. 
The proof can be found in [Q . Now we get 



,,(") 



(g; q)k{a'^q^, abq^+^, a-'^bq-^] q)^ 



j(a) ^ \<^o, q)oo sr^ l-J g ^ " ya q ,a q , gjoo r, (a).L-i /a 7^ 



fc=0 

When b also becomes larger than one, there appear the terms which come from poles at 
Z2 = bq^ and Z2 = {bq^Y^ in (|0|). When (pOI) and 

b>bq>bq^> ■■■> bq""^'^ > 1 > 6g"'"+^ > • • • (A.8) 

hold, we have 

I, = jf ) + /(«) + jf '0) + jf '^) + jf ''^) , (A.9) 

jm ^ iab;q)oo{q;q)l y^ (-)^gM^-i)/^(5V^ r^g-^^ g)^ 

' 26 ^ (g;g)fc(6V;g)oo 

^ f _^ izlzf;qUil + z,)il + zY)Y-' ^^^^^^ 

Jco 27[izi {azi,az^^,bq^+^Zi,bq''+^zY,b-^q-^Zi,b-^q-^zY]q)oo' 

' ~ hmh-1 KT,q)k{q;q)m-k-i 

(g, ab, b'^q^^, fc-^g"^^, fe^g^'", fe-^g-^'"; g)^^ ^ 

^ (g™-'=,62gfc,52g™+fc+1^5-2^-fc-m^gj5^m^^5-l^-m.^)^' >> ' '' 

^(b,a) _ y^ .r- (_)fc+»^gfc(fc+l)/2+n^(m+l)/2[^W]L-l 

r^ .4^ ™ Kq;q)k{q;q)m 



{q,ab,b'^q^^,b ^g '^^,a?q^'^,a ^g ^™;g)c 



(A.12) 



(fo^g/^, a^g"*, abq''+"^+^, a-ifog'^+i"™, ab-^q'^-'', a-ife-ig-'^"™; g)oo ' 
Lastly, when a < 1 and (|A.8| ) holds, we have 

/^ = jf)+/M) + /f=^). (A.13) 

Now we turn to the calculation of the asymptotic expression of the integral Ji. 
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The case a < 1 and b < 1 

In this case, Ji = /} . We evaluate the asymptotic behavior of /} by employing the 
steepest decent method. First we change the variable from 6, ip to u,y as 

l + cose = 2e-"/^, (A.14) 

l + cos(^ = 2ye~"/^, (A.15) 

to obtain 



(„2ie -2ie 2iLp p-2iip. \ 

X ^ ' ' ' '^^°° CA 16) 

(ae*^ ae-'^ ge^(^+v'), ge-*(^+¥'), ge^(''-'^), ge-*(^-'^), fee*"^, fee-*"^; g)oo ' ^ ' ^ 



Here Q and y? are considerd as functions in u and y through (|A.14| ) and ( |A.15| ) 
respectively. Now we can take the limit L — > oo in the integrand. Changing the 
variable y back to </? by 

l + cosv? = 2?/, (A. 17) 

we have 

v/?(g;g)L4^+i /" (e2^^e-2^^;g)oo 



/i ~ - " ^"'"^°° 3 / d<^^ . ^ .' .'":°° , . — . (A.18) 



Using the formula ( [A. 61 ), we get 



Ji ~ ^ ^^. (A.19) 

(a,6;g)^L2 

The case a < 1 and 6 > 1 

The integral Ji is given by ( |A.13| ). The main contributions come from /} and 
/} ' . Each contribution behaves as 4^. whilest the normalization Zi^ behaves as 
[(1 + 6)(1 + 6^"'^)]^ for this case. Since (l + 6)(l + 6~^) is larger than 4, I\ is negligible 
compared to Z^. So we do not compute the explicit expression. 

The case a > 1 

The main contribution comes from the /c = term in the summation of JJ" . We 
have 

a(l - ah){a-\ q)Ji\ + a)(l + g-^)]^'^ 
h ^ 1 _i, s ■ (A.20) 
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Appendix B Evaluation of Integral I2 

In this Appendix, the integral /2, 
/2 = -(a6;g)oo(g;g)c 



4 J ZTTZZi J Z7r^Z2 

(02:1, az^^, qziZ2, az^^z^^, qziz^^ , qz^^z^^, bz2, bz^^; q)^{z2 + z^^ - Zi - z^^) ' 

is evaluated. First, for the case where a,b < 1, both of the contours of zi and Z2 are unit 
circles. We have 

/ =/(o)= r^ r^ ie^'',e-^'',e''''^,e-'''^-q)oo 

^ ^ io 27r Jo 27r (ae*^ae-*^ge*(^+'^),ge-*(^+'^),ge*(^-^),ge-*(^-^),6e*'^,6e-*'^;g)oo 



^ j2(l + cosg)rH2(l + cos^)]^-^- ^ 
2 cos if — 2 cos ^ 

Since we will calculate the average density profile near the right boundary, we set / = 
L — j + 1. Similarly to the case of the integal Ji, when a or 6 or both of them become 
larger than one, there apper other contributions besides I2 ■ When ( |A.3|) and (|A.8| ) hold. 
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the result is 

/2 = 4°) + 4'^'°) + /('^•'^) + /("'^) + jf -o) + jf ^) + jf «), (B.3) 



2a ^ {q\q)k{a^q^\q)c 



2 o„ / > /„. „\ /„9„fc. „\ ' 






/oo 



"*"■' / NA:+m^fc(fc+l)/2+m(m+l)/2r\{a)lL-ir\(6)li-/ 

.(a,fe) _ \^ V^ \~) q '" [Afc J [Am \ 



{q,ab,a?q^^,a ^g "^^ ,b'^q^™-,b ^g ^™;g)c 



(a^gfc, 62gm^ ahq^+'^+^, afe-^g^+i"™, a-^bq"'-'', a~^b~^q-''-"^; g) 



) y/oo 






(_)A.gfc(fc-l)/2(^2^2/c^ 5-2^-2fc. g)^[Af ]'-! 



"^ ' "^ ' / \k^k(k-l)/2+(m-k)(m-k-l)/2\\(b)]l-l\\{b)]L-l 
j(b,b) ^ _ V^ V^ (-J g ^ ^^ ^ ^^ " [Afe J [Am J 

2 ^L^ ^-^ hin-n\,in-n\ , , 



A;=0 m=k+l 



5 

loo 



'r-j: E 



^''^ '_^|fc+m^fe(fe+l)/2+m(m+l)/2rJ)^{b)]/-l|-)^{a)]L-^ 



KM Y- Y- I- J V('^+l)/2+m(m+l)/2[^W]/-l[^ 



{q,ab,b'^q^'',b ^q ^'',a^g^'",a ^g ^™;g)c 



{b^q'', a^q"", abq''+"^+^, a-^bq^+^-'^, ab-^q'^''', a-^fe-ig-^"™; g)oo ' 



(B.6) 



2& Xo {q;q)k{b'^q'';q)oo 

J Co 2vr2Zi {azi,az^'^,bq''+'^zi,bq^+'^z^^,b-^q^^zi,b-^q^^z^'']q)^' 



{q,ab,Pq'^'',b ^q ^'',6^g^™,6 ^g ^"'; g)oo ^^^ c 

^ {q"'-'',b'^q^,b^q"'+''+^,b''^q~''-'^,abq"',ab-^q-"';q)c ' 



(B.9) 



Now we consider the asymptotic expression of the integral /2. We take the limit L ^ oo 
at first and then take the limit / — *> oo. 

• The case a,b < 1 

In this case, I2 = 4 • The evaluation for this case proceeds analogously to the 
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r(0) 



evaluation of /{ . Changing the variables 9, ip to u,y as in (|A.14| ) and ( |A.15| ), we get 



Lf 



duu^e "^ 



1 - e-"/^ 



-M/-L 



dy 






X- 



f'„2J6» p-2i6» „2i(/3 p-2iip. 



g)c 



(B.IO) 



(ae^^ ae-*^ ge^(^+'^), ge-*(^+'^), ge*(''-'^), ge-^(^-'^), fee^-^, fee"^"^; g)oo ' 

We take the limit L ^ cxd in this expression and change the varialbe y back to ip to 
get 



v^(g;?)L4 



2 /iL+2 /"TT 



2(a;g)^L3 Jo 



n ^ 



1 + COSty?) 



(ge^*^, ge-2*^; g) 



(ge*'^, ge*'^, ge"*"^, ge"*'^, 6e*'^, 6e*'^; g)oo ' 

(B.ll) 

We can consider the limit / ^ oo by using the steepest descent method. We have 

2 ■ 4^+i7r 



{a,b]q)l^L3h 



(B.12) 



The case 1 < a < g ^ and 6 < 1 

In this case, the integral I2 is given by I2 = I2 +-^2 ' • The main contribution comes 
from the k = term in L 



{a,0) 



h ^ -{ah, a-2; g)^(g; g)^[(l + a)(l + a-^)]^"' 

dz2 {zl Z2^- g)oo[(l + 2:2) (1 + z^^)]'-^ 



X 



^c-Q 27riz2 (a ^Z2,a ^z^ ,aqz2,aqz2 ,bz2,bz2 ; g)oo 
Taking the limit / ^ 00 by using the steeptest descent method, we get 

4'(a6, a-2; g)^(g; g)^[(l + a)(l + a^^)]^-' 



)2 

/oo 



(B.13) 



(B.14) 



a^/^Th{a~^, qa, b; 

The case a > q^^ and aq > b 

In general, the integral for this case has the expression ( [B.3| ). The main contribution 



comes from the k = 0,m = 1 term in /2 '" . We have 

-2^-2 



^i^-i\i/-i 



1 _ ab)il - a-V')(a-'; g)oo[(l + a)(l + a-^)]^-^[(l + ag)(l + a'^-^)] 



(a ifeg i;g)c 



(B.15) 



The case a > b > aq and 5 > 1 

The ] 

have 



The main contribution for this case comes from the k = 0,m = 1 term in Jg ' • We 



(1 - ab)iq, a-\ b''; g)^[(l + a)(l + a-^)]^-'[(l + 6)(1 + b'')] 
a{a~^b, a~^b~^, ab~^q; g)oo 



-IM/-1 



(B.16) 
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The case b > 1 and b > a 

In this case, the normahzation Z^ behaves as [(1 + 6)(1 + b^^)]^. All contributions 
to I2 can be neglected compared to Z^. Hence we do not compute the asymptotic 
expression explicitly for this case. 
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Figure Captions 



Fig. 1 : One-dimensional partially asymmetric simple exclusion process with open bound- 
aries. Particles have hard-core exclusion interaction and tend to hop to the right (resp. 
left) nearest neiboring site with rate pr (resp. pl)- There are also particle injection (resp. 
ejection) at the left (resp. right) edge. 



Fig. 2 : Space-time diagram of the ASEP from Monte-Carlo simulations. The holizontal 
axis represents the site number j whereas the vertical axis represents time. The existence 
of pariticle is represented as a black point. The lattice length is taken to be L = 200. The 
bulk hopping rates are taken io he pr = 1^ pl = 0. After some transient time, the system 
practically goes to a steady state. The steady state depends crutially on the values of the 
boundary parameters. 



Fig. 3 : Average density profile of the ASEP from Monte-Carlo simulations. The holizontal 
axis represents the site number j whereas the vertical axis represents the average density. 
The lattice length is taken to be L = 200. The bulk hopping rates are taken to be 
Pi? = 1, PL = 



Fig. 4 : The phase diagram of the current. Regions A, B and C are called the low-density 
phase, the high-density phase and the maximal current phase respectively. 



Fig. 5 : The phase diagram of the correlation length in the a-j3 plane for the partially 
asymmetric case. The low-density phase (resp. high-density phase) is divided into three 
phases, Ai,A2 and A^ (resp. -81,-82 and -B3). 



Fig. 6 : The correlation length ^ for the case pn = ^,Pl = 0.9, /3 = 1. The solid line 
is the theoretical 
simulation data. 



is the theoretical prediction given by ^ = 1/ln - J^- ^ whereas the black dots are the 
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3 
{A) low density phase 
PR = 1.0, PL = 0.0 
a = 0.2, (3 = 1.0 



3 
(C) maximal current phase 
PR = 1.0, PL = 0.0 
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coexistence line 
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(B) high density phase 
PR = 1.0, Pi = 0.0 
a = 1.0, /3 = 0.2 
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